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ABSTRACT 

We describe a fast and accurate method for estimation of the cosmic microwave background 
(CMB) anisotropy angular power spectrum — Monte Carlo Apodiscd Spherical Transform 
EstimatoR. Originally devised for use in the interpretation of the Boomerang experimental data, 
MASTER is both a computationally efficient method suitable for use with the currently available 
CMB data sets (already large in size, despite covering small fractions of the sky, and affected by 
inhomogeneous and correlated noise) , and a very promising application for the analysis of very large 
future CMB satellite mission products. 



1. Introduction 

During the past decade since the ground-breaking 
discovery of the cosmic microwave background ra- 
diation anisotropy by the COBE satellite (Smoot 
et al. 1992), numerous successful measurements of 
microwave sky structures have provided us with the 
data for powerful tests of the current cosmological 
paradigm, and created an unprecedented opportu- 
nity to estimate key parameters of the candidate 
theoretical models of the Universe. 

Recent ground-based and balloon-borne experi- 
ments with improved sky coverage, angular resolu- 
tion, and noise performance (see de Bernardis et al, 
2000, Hanany et al, 2000, Padin et al, 2000, Jaffc 
et al, 2001, Lee et al, 2001, Halverson et al, 2001, 
Pryke et al, 2001 and references therein for some 
of the most recent experiments and their interpre- 
tation) have both given us a taste of what future 
satellite missions MAP 1 and Planck 2 should accom- 
plish, and revealed the growing challenges that we 
will have to meet in the analysis of the forthcoming 
CMB data sets. 

In the currently favoured structure formation 
model of inflation induced, Gaussian distributed, 



Microwave Anisotropy Probe, http://map.gsfc.nasa.gov/ 
http: / /astro. estec.esa.nl/Planck/ 



curvature perturbations all the statistical informa- 
tion contained in a CMB map can be summarised 
in its angular power spectrum Ci. General maxi- 
mum likelihood methods for extracting Ci from a 
Apix-pixel map with non uniform coverage and cor- 
related noise (Gorski 1994, Bond 1995, Tegmark 
& Bunn 1995, Gorski 1997, Bond, Jaffe, & Knox 
1998, Borrill 1999b) involve computations of com- 
plexity ~ -/Vpix 3 , and become prohibitively CPU ex- 
pensive for the -/Vpi x > 10 4 maps produced by cur- 
rent experiments. With the presently anticipated 
computer performance such methods appear totally 
impractical for application to the -/Vpi x > 10 6 maps 
expected from the future space missions (Borrill 
1999a). Hence, there is a well recognised need for 
faster, more economical, and accurate Cg extraction 
methods, which should enable a correct comological 
interpretion of the CMB anisotropy observations. 

In this paper we introduce and discuss a new 
method for fast estimation of the CMB anisotropy 
angular power spectrum from fluctuations observed 
on a limited area of the sky. This method is based on 
a direct spherical harmonic transform (SHT) of the 
available map and allows one to incorporate a de- 
scription of the particular properties of a given CMB 
experiment, including the survey geometry, scan- 
ning strategy, instrumental noise behaviour, and 



possible non-gaussian and/or non-stationary events 
which can occur during the data acquisition. The 
estimated power spectrum is affected by the un- 
wanted contribution of the instrumental noise and 
the effects of any necessary alteration of either the 
recorded data stream (such as high pass filtering) or 
the raw map of the observed region of the sky, which 
are introduced during the data analysis. These ef- 
fects are calibrated in Monte Carlo (MC) simula- 
tions of the modeled observation and analysis stage 
of the experiment and can then be removed, or cor- 
rected for in the estimated power spectrum. The 
harmonic mode-mode coupling induced by the in- 
complete sky coverage is described analytically by 
the SHT of the sky window and corrected for in or- 
der to obtain an unbiased estimate of the Cg. Here- 
after we refer to this method with an acronym MAS- 
TER ( Monte Carlo Apodised Spherical Transform 
EstimatoR) . 

Netterficld et al. (2001) described an application 
of this method in the extraction of the CMB an- 
gular power spectrum Cg, for 75 < I < 1025, from 
the sky map (analysed region comprised ~1.8% of 
the sky covered with 57000 pixels of 7' size) made by 
coadding four frequency channel data of the 1998/99 
Antarctic long duration flight of the Boomerang 
experiment (Boom-LDB). The first derivation of 
the CMB anisotropy spectrum from the same data 
(de Bcrnardis et al. 2000) involved the MADCAP 
method (Borrill 1999b) applied to a smaller sub- 
set of the data (one frequency channel, ~1% of the 
sky covered with ~ 8000 pixels of 14' size). The 
MADCAP approach is too CPU intensive for re- 
peated applications to the new, enlarged subset of 
the Boomerang data, and, hence, the MASTER 
approach was the method of choice for extraction 
of the h\gh-£ angular power spectrum of the CMB 
anisotropy. 

Other fast methods have recently been proposed 
for estimation of the angular spectrum of the CMB 
anistropy. Szapudi et al. (2001) advocate the use 
of the 2-point correlation function for extraction of 
the angular power spectrum from the CMB maps. 
The computational demands of this method scale 
quadratically, <~ iV p j x 2 , with the size of data set 
(that may be improved to ~ N p i x log N p ^) . In 
the same way as in the case of MASTER, the ef- 
fects of the noise and correlations of the derived 
Cg-s are quantified by Monte Carlo simulations (al- 
though the demonstrated applications involved only 
the case of a uniform white noise). Dore, Knox & 
Peel (2001b) proposed a hierarchical implementa- 
tion of the usual quadratic Cg estimator with a com- 
putational scaling proportional to A p ; x 2 , that may 
be reduced to ~ A p ; x (with a large prefactor) at the 
price of additional approximations. 



Experiment specific techniques have also been 
proposed : Oh, Spergel, and Hinshaw (1999) de- 
scribed a fast power spectrum extraction technique 
designed for usage with the MAP satellite data. 
Their method scales like ~ A p ; x 2 with the size of 
the pixellised map, and takes advantage of uncor- 
rected pixel noise with approximate axisymmetric 
distribution on the sky. Wandelt (2000) advocates 
the use of the set of rings as a compressed form of the 
Planck data set from which to extract optimally the 
Cg-s in the presence of correlated noise. The appli- 
cability of this approach is limited by its assumption 
of the the symmetry of the scanning strategy. 

This paper is organised as follows: In section 2 
we describe how a data stream of observations is re- 
duced to a CMB fluctuation map, and how the an- 
gular pseudo power spectrum Cg is extracted from 
such a map by SHT. In section 3 we show how an un- 
biased estimate of the true underlying power spec- 
trum can be recovered from the Cg-s with the aid 
of the Monte Carlo simulations. The tests of the 
method on simulated Boom-LDB observations are 
described in section 4, and the application of the 
method is discussed in section 5. 

2. Prom Time Ordered Data to Pseudo 
Power Spectrum 

Single dish CMB experiments produce for each 
detector a data stream, or the time ordered data 
(TOD), of the direction of observation and the sky 
temperature as measured through the instrumental 
beam. We assume that the beam is known, that it 
is close to isotropic in the main lobe, that the side 
lobes are negligible, and and that the pointing at 
each time is known to an accuracy better than the 
size of the main lobe of the beam. Exceptions to 
these assumptions will be addressed in section 3.6. 
We will also assume that all the TOD samples af- 
fected by transient events, such as cosmic ray hits, 
have been removed and that in order to preserve the 
TOD continuity the resulting gaps are filled with 
fake data having the same statistical properties as 
the genuine observations (eg, Prunet et al, 2000, 
Stompor et al, 2000). 

2.1. From TOD to Sky Map 

The data produced by each detector at a time t 
can be modeled as 

dt = PtpAp + n t , (1) 

where A p is the sky temperature, that we assume 
to be pixeliscd and smoothed with the instrument 
beam, P tp is the pointing matrix, p is the pixel index 
and n t is the instrumental noise. 



If the TOD noise is Gaussian distributed with a 
known correlation function N t t> = (n(t)n(t')), the 
optimal solution for the sky map 

m p = {Pl t N-}P t , p ,)^Pl t N-}d t , (2) 

minimises the residual noise in the pixellised map, 
A p - m p (Lupton 1993, Wright 1996, Tegmark 
1997). While being completely general this pro- 
cedure is impractical for very long TOD streams 
because of the required inversion of the large ma- 
trix N tt i. A simplification is possible under the as- 
sumption of the TOD noise being piece-wise sta- 
tionary, and its correlation matrix representable as 
circulant, N tt > — N(t — t 1 ). Eq. (2) can then be 
solved either directly (with a computational scaling 
of ~ 7Vpi x 3 ), or by using iterative methods as dis- 
cussed by Wright (1996), or Natoli et al. (2001) (in 
which case the computation time is dominated by 
Fourier space convolutions of the TOD correspond- 
ing to the product N^, dp in Eq. 2). Iterative ap- 
proaches scale like Ni ter N T log N T , where N T is the 
number of time samples, and Ai tor is the number of 
iterations. -/Vj ter depends on the required accuracy 
of the final map, and it is of the order of a few tens 
in the case of a conjugate gradient method of linear 
system solution (Natoli et al. 2001). 

If the TOD noise properties are not known be- 
forehand, however, as is generally the case, the 
Eqs (1) and (2) can be solved iteratively together. 
This returns at each time step an estimate of the 
noise stream, n(t), and, hence, of the noise time 
power spectrum. The required computational scal- 
ing involves a somewhat larger Mter (Ferreira and 
Jaffe 2000, Prunct et al. 2000, Stompor et al. 2000, 
and Dore et al. 2001a). 

Since the MASTER method requires repetitive 
TOD simulations, processing, and map making, the 
iterative solution of Eq. (2) can be too time consum- 
ing for practical applications. Therefore, to avoid 
the necessity to iterate, we use a suboptimal, fast 
map making method involving the high pass filter- 
ing of the TOD stream, which improves the long 
time scale behaviour of the noise, and reduces the 
striping of the resulting map (see section 3.2). The 
map solution is now 

m p = (NoUp))- 1 E P Ptf(t - 0*'. (3) 
w 

where N b s (p) = P^ t P tp is the number of observa- 
tions in the pixel p, and / denotes the high pass 
filter. The computational scaling is now reduced 
to iV r logAr T . Clearly, Eq. (3) is only equivalent 
to Eq. (2) if the TOD noise is white, i.e. N a > = 
No5vi Yac (t — t') (in which case the filter would be re- 
duced to / — 5r>i iac (t), i.e. no filtering would be ap- 
plied). While the application of the high pass filter 



reduces the long term noise correlations, it degrades 
the CMB signal at low frequencies (see Fig. 1) and 
affects the resulting angular power spectrum derived 
from the filtered map solution Eq. (3). This effect is 
quantified and corrected for with the Monte Carlo 
simulations and analysis involving the filtered map 
making technique applied to the simulated TODs of 
the pure CMB signal. This procedure will be dis- 
cussed in detail later on. 

2.2. From Sky Map to Pseudo Power Spec- 
trum 

A scalar field AT(n) defined over the full sky can 
be decomposed in spherical harmonic coefficients 

a tm = J dnAT(n)r; m (n), (4) 

with 

e 

AT(n)=^^a <m y <ra (n). (5) 

£>0 m=-l 

If the CMB temperature fluctuation AT is as- 
sumed to be Gaussian distributed, each ai m is an 
independent Gaussian deviate with 

{aim) = 0, (6) 

and 

{ai m a* t , m i) = 6u'5 mm >{Ct), (7) 

where {Ce) = Cf 1 is specified by the theory of 
primordial perturbations, and parametrised accord- 
ingly, and 5 is the Kronecker symbol. An unbiased 
estimator of Cj h is given by 

1 1 

° e = 21+1 ^ \ aim \ 2 ■ ( 8 ) 

m= — I 

C*£-s are x^-distributed with the mean equal to C| h , 
v = 2£ + 1 degrees of freedom (dof), and a variance 
of 2Cf 2 jv. 

In the case of CMB measurements the tempera- 
ture fluctuations can not be measured over the full 
sky, either because of ground obscuration or galac- 
tic contamination for example, and a position de- 
pendent weighting W(n) can also be applied to the 
measured data, for instance to reduce the edge ef- 
fects. If / s ky represent the sky fraction over which 
the weighting applied is non zero, then 

UyWi = J- / dnW*(n) (9) 

471 " J47T 

is the i-th moment of the arbitrary weighting 
scheme. The window function can also be ex- 
panded in spherical harmonics with the coefficients 



w £m = J dnW (n)Y f * m (n) , and with a power spec- 
trum 

^ = ^TiEi w ^i 2 < ( 10 ) 

m 

for which W{£ = 0) = 47r/ sky 2 wf and 
^ e > W(£)(2£ +1)=4tt f sky w 2 . 

A sky temperature fluctuation map AT(n) on 
which a window W(n) is applied can be decomposed 
in spherical harmonics coefficients 

a im = J dnAT(n)W(n)y/ m (n) (11) 

« fi p ^AT(p)W(p)F/ m (p), (12) 
p 

where the integral over the sky is approximated by 
a discrete sum over the pixels that make the map, 
with an individual surface area fL. 



The pseudo power spectrum C can be defined as 
1 



21+1 E ' 

m=— I 



(13) 



3. Prom Pseudo Power Spectrum to Full 
Sky Power Spectrum Estimator 

The pseudo power spectrum Ce rendered by the 
direct spherical harmonics transform of a partial sky 
map, Eq. (13), is clearly different from the full sky 
angular spectrum Ce, but their ensemble averages 
can be related by 



(C t ) = ^2Mu'{Ct>), 



(14) 



where M«> describes the mode-mode coupling re- 
sulting from the cut sky. As described in the ap- 
pendix, this kernel depends only on the geometry of 
the cut sky and can be expressed simply in terms of 
the power spectrum We of the spatial window ap- 
plied to the survey (see Eq. (A31) and Eq. (A14) for 
the spherical and planar geometry, respectively). 

The effect of the instrumental beam, experimen- 
tal noise, and filtering of the TOD stream can be 
included as follows 

(Ce) = Y, M u >F fJ Bl(C t ) + (N t ), (15) 



The computation of Eq. (12) for each (£,m) up 
to £ = £ max performed on an arbitrary pixelisation 
of the sphere would scale as A^ p i x £ max 2 . However, 
if one uses an adequate lay out of the pixels to ex- 
ploit the symmetries of Spherical Harmonics, such 
as for example the ECP (Muciaccia et al. 1997), 
HEALPix (Gorski et al. 1998), or Igloo (Crittenden 
and Turok 1998) this computation actually scales 
like A p i x 1 / 2 ^ max 2 . In our implementation of the 
MASTER method, after application of the window 
function on the map, the program anafast from the 
package HEALPix was used to compute the pseudo 
power spectrum. 

Wandelt, Hivon & Gorski (2000) showed that the 
marginalised likelihood P(Ce\C th , N) of the pseudo 
Ce for a given underlying theory C th and a given 
noise covariance N can be computed analytically 
in C(A p i x 1 / 2 £ max 2 ) operations, under conditions of 
axisymmctric sky window function and (non nec- 
essarily uniform) white noise. Under these assump- 
tions P(Ce\C th ) can be used to perform a maximum 
likelihood fit of the cosmological parameters to the 
observed data set. Hansen et al (2001) extend this 
approach by using the full pseudo Ce covariance ma- 
trix. We will now build, starting from the measured 
Ce, and under more general conditions on the noise 
properties and shape of the observing window, a 
new estimator of the full sky power spectrum that 
can be compared directly to C th . 



where Be is a window function describing the com- 
bined smoothing effects of the beam and finite pixel 
size, (Ne) is the average noise power spectrum, and 
Fe is a transfer function which models the effect of 
the filtering applied to the data stream or to the 
maps. The determination of each of these terms 
will be described below. 

It is often assumed that the xt=2t+i distribu- 
tion of Ce-s on the full sky can be generalised to 
cut sky observations by scaling v to the number of 
dof effectively available. Given the large value of v 
the central limit theorem is also invoked to further 
simplify this to a Gaussian of the same mean and 
variance. From these successive (and excessive) sim- 
plifications we will only retain, as a rule of thumb, 
that the rms of Ce averaged over a range A£ is ap- 
proximately 



AC, 



Ce+ BW)) 



(16) 



with 



ue = (2£ + l)A£f sky ^, (17) 

where the factor w\jw^ accounts for the loss of 
modes induced by the pixel weighting. We will 
show in section 4 how this compares to the results 
of Monte Carlo simulations. 

3.1. Mode-mode Coupling Kernel 



The resolution in £ of the measured power spec- 
trum is ultimately determined by the extent of the 




Fig. 1. — Simulation of the Boom-LDB experiment and application of MASTER to extract the CMB angular 
power spectrum. The oval contour on the maps shows the ellipse (distorted by projection) within which the 
power spectrum is estimated (/ s ky = 1.8% of the sky). Top left panel: A random realisation of the CMB sky 
from the theoretical model described by the power spectrum shown in the top right panel (red line) . Middle left 
panel: A noiseless map of the same region of the sky made from the TOD with actual Boom-LDB pointing and 
processed with the 100 mHz high-pass filter (see text). Bottom left panel: the difference between the two CMB 
sky maps shown above, which shows the component of the sky signal lost due to the combination of Boomerang 
scanning and data processing. Middle right panel: A simulation of the same Boomerang CMB sky map with 
the instrumental noise included. Bottom right panel: Integration time per pixel for the actual scanning of the 
Boom-LDB channel B150A; the average integration time is about 500s/Deg 2 . Top right panel: The input power 
spectrum smoothed by the beam and pixel window function (red line); The average angular power spectrum 
of the instrumental noise (black line); The pseudo Cg-s directly measured on the sky map shown in the middle 
right panel and divided by / s k y (orange line); The binned MASTER estimate of the full sky power spectrum 
after removal of noise contribution and correction of the effect of the high pass filtering and mode- mode coupling 
(blue histogram). 



oberved area of the sky, its geometrical shape, and 
the pixel weighting W(n) applied to the survey (see 
Hobson & Magueijo 1996, and Tegmark 1996). Al- 
though we only tested numerically the method on a 
circular or elliptically shaped window, nothing pre- 
vents the use of a more complex window, specially 
for a pixel starved experiment with a nonconvex 
survey area, for which a well designed apodisation 
could help improving the achievable spectral reso- 
lution. We will show in section (4) how the choice 
of window changes the estimated Cg spectrum. 

3.2. TOD Filter Transfer Function 

The transfer function Fi introduced in Eq. (15) 
describes the effect of any filtering applied to the 
TOD stream or to the map on the angular power 
spectrum. A specific example of the latter is the 
removal of parallel stripes extending along a direc- 
tion different from the scanning direction observed 
in some channels of the Boom-LDB data (Netter- 
field et al. 2001, Contaldi et al. 2001). The filtering 
of the TOD has broader applications, however, and 
can take the form of a high pass filtering that serves 
several purposes, as follows: 

• reduce the contribution of the low frequency 
noise (1/f noise) to the map, specially if the 
scanning strategy and/or the map making 
technique used do not optimize the removal 
of these modes, 

• reduce the scan or spin synchronous noise, 
which may appear at the scan frequency and 
its harmonics, 

• remove from the signal the contribution from 
the large scale anisotropies, which are poorly 
constrained on an incomplete sky survey, and 
are likely to contaminate the estimated power 
spectrum at all the smaller scales. 

It should be noted that the validity of Eq. (14) 
for any sky window relies on the fact that the sta- 
tistical properties of the full sky temperature fluc- 
tuations are isotropic, as implied by Eq. (7). This 
assumption is broken by the high pass filtering of 
the data stream which creates a preferred direction 
on the sky. Therefore the introduction of the single 
scalar function Ft in Eq. (15) should be seen as a 
simplifying ansatz for a more complex reality. We 
can numerically test the validity of this ansatz by 
showing for instance that Fg depends weakly on the 
choice of the underlying power spectrum. 

Given an input CMB power spectrum C\ , a 
number, of noise- free full sky realizations of 

this spectrum can be simulated (using for example 
the program synfast from HEALPix). These maps 
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Fig. 2. — Transfer function of the measured angu- 
lar power spectrum corresponding to the sharp high 
pass filtering at 100 mHz of the Boom-LDB TOD is 
shown here derived for for two different input CMB 
anisotropy power spectra. The upper panel shows 
the two input spectra corresponding to a flat Uni- 
verse (green solid line), and an open Universe (red 
dashed line). The bottom panel shows the respec- 
tive transfer functions, which differ by less than 1% 
up to £ = 1000. The black dotted line shows the an- 
alytical prediction for a toy model of parallel scans 
(Eqs. B14, B15). 

are then "observed" using the actual scanning strat- 
egy, the resulting model TODs projected back on 
the sky using Eq. (3), and their individual power 
spectra Cg are extracted using Eqs. (12) and (13). 

Equation (15) can then be applied to the aver- 
age (Ci)mc of these measurements to determine 
the transfer function. In order to avoid inverting 
the kernel Mw this system can be solved itera- 
tively. In the case of a scanning experiment, such 
as Boomerang, Appendix B shows how a model of 
parallel scan can render a first order solution F^ . 
In the more general cases, a starting solution could 
be F e {0) = (C i )Mc/(UyW 2 B$Cf i y If (Ci) M c is re- 
placed by its running average S((Cg) mc) (typically 
computed over A£ = 50 points) it appears that one 
iteration is sufficient to obtain a stable estimate of 
the transfer function 

F e = Ff ] (18) 
S ({C t )Mc) - Ee M ltl F^B%{C e ) 

We used Eq. (19) to compute the transfer func- 
tion in Monte Carlo simulations of Boom-LDB ob- 



servations on an clliptically shaped region which 
comprises 1.8% of the sky (see section 4 for details) 
for two different input power spectra; the first one 
corresponding to a flat universe with a first peak 
at ^p C ak — 220, and the other one corresponding 
to an open universe with £ pca k — 420. Figure 2 
shows that the resulting Fg-s are almost identical, 
with a discrepancy smaller than 1% in the range 
50 < £ < 1000. This justifies the use of the sim- 
ple ansatz (15) as a model of the effect of the time 
filtering on the angular power spectrum and demon- 
strates that the determination of the transfer func- 
tion can be done nearly independently of any as- 
sumptions about the actual CMB power spectrum. 

Using approximation (16) one expects that the 
error 8Fg done on the MC estimation of the transfer 
function decreases as 




1 



(2^+l)A£/ sky 



(b) 



(19) 



MC 



So if A^/ sky ~ 1 and £ > 100, an estimate of the 
transfer function better than 1% can be obtained in 
4'c £ 100 realisations. 

The computation of the transfer function Fg is 
required because of the filtered map making tech- 
nique used (Eq. 3), which alters the signal in low 
frequency modes. However, even if a more sophisti- 
cated map making were used, the map obtained is 
usually not an unbiased representation of the true 
sky because of various systematic effects present in 
the data, and the computation of Fg is still neces- 
sary. 

If several detectors with different beams are anal- 
ysed simultaneously, the realisation of the same sky 
with a different smearing can be used as an input to 
the TOD simulation of each detector. In the anal- 
ysis of the coadded map an effective beam window 
function B(£) can be modeled as the weighted av- 
erage of all individual beams (see Wu et al. 2001). 
If the actual effective beam were different from this 
model the difference would be reflected in the trans- 
fer function F. 

3.3. Noise Power Spectrum 

From the point of view of modeling the CMB 
experiment the simplest form of instrumental noise 
is the stationary, white, Gaussian process. Real- 
ity however is usually not as simple, the actual ex- 
perimental noise is often non-stationary, "coloured" , 
sometimes non-Gaussian, and correlated to some in- 
ternal variables of the instrument, such as its accel- 
eration, the cold plate temperature, the orientation 
relative to the sun, to the balloon or to the ground. 
For many of those reasons the noise often can not 
be efficiently averaged out. However, if these noises 



can be modeled to a reasonable accuracy, they can 
be included, at little or no extra cost, in the Monte- 
Carlo pipeline described here, and their effect on 
the measured Cg-s can be assessed, and possibly re- 
moved. 

As mentioned in §2.1 an estimate of the noise 
time correlation function N tt i and its time power 
spectrum can be extracted from the actual data 
stream. Using this information a fake Gaussian 
"noise stream" can be simulated and projected on 
the sky with the actual scanning strategy using 
Eq. (3), with the same high pass filtering /. The 
power spectrum N{£) of the resulting noise map is 
extracted according to Eqs. (12), and (13) and the 
whole process is reproduced as many times as nec- 
essary to obtain a Monte-Carlo estimate of the av- 
erage noise angular power spectrum, (N(£))mc- If 
several detectors are analysed simulatenously, and 
their noise is known to be correlated, these corre- 
lations can be included in the noise stream simula- 
tions. 

3.4. Estimated Power Spectrum 

In order to reduce the correlations of the Cg-s in- 
duced by the cut sky, and also to reduce the errors 
on the resulting power spectrum estimator, it is con- 
venient to bin the power spectrum in £. The slowly 
varying "flattened" spectrum Cg = £(£ + l)Cg/2n is 
a preferable candidate for such binning. For a set of 
fibins bins, indexed by 6, with respective boundaries 

4ot < 4dgh < Aow^' onc can define the binning 
operator P as follows 



Pu 



2tt _ /<0 ' 

'-low How 

= 0, otherwise 



if 2 < 



t {b) 

low 



<£< 



f(b+l) 
low 

(20) 



and the binned power spectrum is Cb = PuCg. The 
reciprocal operator (corresponding to a piece-wise 
interpolation) then reads 



Qtb = 



2tt 



if 



= 0, 



otherwise. 



2 < e < 



< 



f(b+l) 
■-low 



(21) 



The two operators above are defined for the flat 
band, disjoint bins. They can be easily modified to 
account for an £ dependent weighting within each 
bin (for example designed to enhance the less noisy 
multipoles) without changing the rest of the discus- 
sion. 

Rewriting Eq. (15) as 

(Cg) = Kgg,(Cg,) + (Ng) (22) 
we then look for a solution to 

PbtKw (CV) - Pu ((C t ) - {N t )) ■ (23) 



This system has £ max unknowns for Ubins equa- 
tions. We seek the solution such that {Ct) — t{l + 
l)(Ce)/2ir is a piece-wise constant. If we replace 
(C t ) by Qt b P w (Cf) = Q eb (C b ) then 

(C b ) = K b - b }P b , e ({C e ) - (N t )) , (24) 

where 

Kw = PbtKwQt'b', 

= P M M u ,F e ,B 2 e ,Q e , b ,. (25) 

An unbiased estimator C b of the whole sky power 
spectrum is then given by 

C b = K bb )P V t (p t - (N e ) MC ) . (26) 

An estimator of the noise "on the sky" can also be 
introduced 

N b = K bb }P b , e (N e ) M c. (27) 

Figure 3.4 shows the kernel K bb > and its inverse 
for the configuration described in section 4. 

It should be noted that the binning of the I space 
is performed at the last stage of the analysis and can 
be chosen after the MC simulations are done. 

3.5. Covariance Matrix of the Estimated 
Power Spectrum 

In order to be able to extract the cosmological 
parameters from the power spectrum estimated as 
described above, one needs to know the errors on 
each C b , and the correlations between the bins. This 
information is contained in the covariance matrix of 
the estimated power spectrum, which can be esti- 
mated as follows. A smooth interpolation of C b is 
used as the underlying CMB power spectrum, and a 
new set of Monte Carlo simulations, including both 
signal and noise, as well as all the experiment pe- 
culiarities, is generated, analysed, and reduced in 
the same way as the real data. This generates a 
set of binned power spectrum estimators {C b }. The 
elements of the correlation matrix are defined as 

C bb ' = ((C b - (C;,)mC*) (Cb 1 - (Cb')MC^ 

MC (28) 

The error bars on C b are then given by the square 
root of the diagonal elements of C 

AC b = Cli 2 . (29) 

3.6. The Algorithm 

Assuming that for a given CMB experiment the 
instrumental beam and the time domain power 



spectrum of the instrumental noise are known, the 
process of estimation of the full sky power spectrum 
from the noisy observations of CMB temperature 
fluctuations can be summarised as follows. 

The required tools are : 

(a) a simulation facility for generation of the ran- 
dom realisations of the CMB sky (e.g. synfast 
from HEALPix), 

(bl) a software model of the experiment which sim- 
ulates observations of the sky using the ap- 
propriate scanning strategy (Eq. 1), and gen- 
erates the model CMB signal TOD streams, 

(b2) a noise simulator that can generate random 
realisations of the noise with an appropriate 
power spectrum; possible non-gaussian noise 
features or cross-correlations between detec- 
tors should be added at this stage, 

(b3) a fast map making facility, which implements 
Eq. (3), and accounts for any alterations of the 
observed TOD stream and/or the produced 
map, 

(c) a software to compute the pseudo power spec- 
trum (Eqs. 12, and 13) from a given apodised 
cut sky map (e.g. anafast from HEALPix). 

After the choice of the sky window apodisation 
function is made the procedure of estimation of the 
power spectrum involves the following steps: 

1. Eq. (A31) is used to evaluate the Ct coupling 
kernel Mu>, which accounts for the effects of 
limited sky coverage and apodisation; 

(s) 

2. a number N^c of noise free Monte Carlo sim- 
ulations of the observed TOD (produced with 
(a) and (bl), projected on the sky with (b3), 
and analysed with (c)) are used to estimate 
the transfer function Ft of any filtering that is 
applied to the actual TOD stream; 

3. a number N^q of pure noise Monte Carlo sim- 
ulations of the TODs (made with (bl) and 
(b2), and then projected and analysed with 
(b3) and (c)) are used to estimate the angular 
power spectrum (N)mc of the noise projected 
on the sky; 

4. the experimental TOD is coverted into a map 
using (b3) and its pseudo power spectrum Ct 
is obtained with (c); 

5. next a set of l-hvos is defined, and an estimate 
of the underlying full sky binned power spec- 
trum is computed using Eq. (26); according 





Fig. 3. — Binned power spectrum coupling kernel Kbb> and its inverse (absolute values shown, with green 
color indicating the negative elements) for an elliptically shaped top hat sky window covering 1.8% of the 
sky. Binwidth is A£ = 50, except for the last bin, for which A£ = 150. The diagonal elements and the 
£b 1 = 200, and 700 rows of both matrices are shown in the bottom panels. 



to Eq. (24) this is an unbiased estimator - 
(Cb) = (Cb)', this will be demonstrated in sec- 
tion 4 with simulated Boomerang observations 
in which the input spectrum is known; 

6. the covariance matrix Cbw (Eq. 28) is com- 
puted from N^, 1 ^ simulations of the whole 
experiment, and the error bars on the binned 
power spectrum are obtained from its diagonal 
elements (Eq.29). 

We assumed the physical beam to be close 
enough to axisymmetric so that its smoothing ef- 



fect on the temperature map is independent of the 
payload attitude along the line of sight. It it were 
not the direct integration of the tempera- 

ture over the beam would have to be performed for 
each time sample in Eq. (1). This operation could 
be very intensive for extremely structured beams, 
unless some symmetries in the scanning strategy, 
such as the one expected for satellite missions, al- 
low for a fast convolution implementation (Wandelt 
& Gorski 2000, Challinor et al, 2000) 

On the other hand, the effect of a pointing inacur- 
racy, whether it is axisymmetric or not, can easily 
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Fig. 4. — Spectral densities of the noise measured 
for one of the Boom-LDB channels is shown for the 
two scanning speeds (1 and 2 degrees per second) 
of the actual CMB observations. These spectra 
were used to generate the simulated observations. 
A sharp high-pass filter was applied to the mea- 
sured and simulated TODs at 100 and 200 mHz, 
respectively. 

be included in the method by modifying Eq. (1) ac- 
cordingly. 

3.7. Computational Scaling 

The algorithm is maximally parallelisable as each 
Monte Carlo cycle of the algorithm can be per- 
formed on a separate CPU. Therefore the total time 
required for completion of the estimation of the 
power spectrum from the single detector CMB ob- 
servations, with -/Vmc = ^mc+^mc+^mc 11 ' cycles 
run on ncpu processors, is given by 

T total = — — -TmC cycle (30) 
"CPU 

with T M c cycle = T( a ) + T( b ) + T (c) , where (a),(b), 
and (c) refer to the CPU time consumption by the 
simulation tools described in section 3.6 (CMB map 
synthesis, observation simulation and map making, 
and map analysis, respectively). In the case of 
joint multi-detector analysis the stage (b) has to 
be repeated for each detector, whereas the stage (a) 
only has to be repeated for each different beam. In 
our own implementation (see section 4 for detailed 
specifications) the MASTER method is executed on 
personal computers (PCs) equipped with 850 MHz 
AMD Athlon CPUs, the following performance is 
achieved (on each processor): 

r w = r (t) , 3 oo s (!^) I,2 G=) 2 , (31) 



and 

where N p i x is the number of pixels over the whole 
sky at the chosen map resolution (not in the cut sky 
map), N T is the number of time samples in the TOD 
set used, and Afft is the typical number of time 
samples on which the required FFTs are computed. 
This leads to an overall CPU time requirement of 

T „„.„, day (^)(_A_). (33 ) 

These timing estimates could be easily improved 
if some repetitive tasks, such as the translation from 
sky coordinates to pixel index in (b), or the evalu- 
ation of Legendre polynomials in (a) and (c) are 
precomputed and stored on disc. 

NB: The most daunting challenge currently fore- 
seen for CMB anisotropy analysis arises in the pro- 
cess of reduction and interpretation of the data that 
will be collected during the ESA Planck mission 
(currently scheduled for launch in 2007). Let us re- 
compute our CPU time requirements to match the 
parameters which describe some of the Planck High 
Frequency Instrument specifications. One chan- 
nel of the HFI with a beam resolution of 5 ar- 
cmin should, during one year of observations, return 
N T ~ 5.6 x 10 9 TOD points, which will be used to 
make a HEALPix full sky map with Apj x = 5 x 10 7 
pixels (of average size of 1.7 arcmin) and which in 
turn should allow us to estimate the angular power 
spectrum of the sky signals up to £ max — 3000. With 
these specifications, and assuming Afft — 1-5 x 10 7 
(one day), the CPU time required for execution of 
our current implementation of MASTER becomes 

Md -(w)(sb) (34) 

on currently available PCs. Since dedicated com- 
puter servers with 32, or more, processors at least 
twice as fast as the PCs that we used are already 
presently available, the quoted total execution time 
for the MASTER method can be reduced to about 3 
days. The extra speed up of the CPUs between now 
and the launch of Planck (a factor of ~ 16 accord- 
ing to Moore's law) should enable the simultaneous 
analysis of several Planck channels, with more so- 
phisticated map making and a larger number of MC 
cycles to improve the power spectrum estimation ac- 
curacy, in a total time of a few days. This means 
that the MASTER approach is fully practical from 
the point of view of demands related to scientific 
analysis of the Planck data. 
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Fig. 5. — Results of the test of accuracy of our method of estimation of the CMB power spectrum are shown 
in application to the simulated Boom-LDB data (with known input power spectrum) with three different sky 
window apodisations (top hat, cosine and Gaussian from top to bottom panels). We used an ensemble of 1352 
Monte Carlo simulations. Left panels: The input power spectrum is shown as a solid black line, and the black 
histogram shows its bin-averaged values. The dashed lines show the average spectra of the simulated noise. 
The red histogram shows our MASTER ensemble mean values of the estimated bin-averaged power spectra, 
which are in perfect agreement (to an error on the mean) with the input bin-averaged theory. For each spectral 
bin three error bars are shown, all centered vertically on the mean values of the MASTER estimates, (Cb/MC, 
and, for clarity, spread around the ^-centers of the bins. The central (red) error bars show the rms values of 
the MASTER ensemble of estimates, while the right-, and left-shifted (green and blue, respectively) ones show 
the theoretical error estimates based on Eqs. (35, 36) derived, respectively, without and with an inclusion of 
the effect of the TOD high-pass filter related spectral transfer function Fi. Right panels: absolute value of 

1/2 

normalized binned spectrum correlation matrix Cw I {CbbCb'b') (see text). Elements with value in [—0.2,0] 
are represented in green, those with value in [0, 0.2] in red and those larger than 0.2 in black. 



4. MASTER Tests on Simulated Boomerang 
Observations 

4.1. Simulations of CMB Observations Dur- 
ing The Long Duration Baloon Flight 
of Boomerang Experiment 

The MASTER method was tested for applica- 
tion to the extraction of the CMB anisotropy power 
spectrum from the data collected by the Boom-LDB 
experiment. These CMB observations comprise a 
total of about 50 10 6 samples which cover about 
4.4% of the sky, and were acquired at two different 
azimuthal sky scan rates of 1 and 2 deg per second 
(for more details on the Boom-LDB flight see Crill, 
2001). 

We have modeled these CMB observations of the 
simulated CMB sky with a scanning pattern iden- 
tical to that of one of Boom-LDB channels. The 
instrumental noise generated in this channel was 
assumed to be Gaussian with a time power spec- 
trum identical to the one measured during the ob- 
servations. The characteristic features of this noise 
power spectrum include a ~ 1// behavior at low fre- 
quency, a knee-frequency of about lOOmHz, a white 
noise level of 130 /iK.s 1 / 2 , a series of lines located 
at the harmonics of the scanning frequency (8 mHz 
at ldps, and 16 mHz at 2dps), and some micro- 
phonic artifacts at 8 Hz (see fig 4). Angular res- 
olution of this channel was FWHM w 10 arcmin, 
and we used the actual measured beam profile in 
the calculations. 

The high-pass filtering applied during the process 
of map-making (Eq.3) was a sharp cut off at 100 
mHz for the ldps scan rate, and 200 mHz for the 
2dps scan rate. 

The sky maps where pixelised using HEALPix 
with 7 arcmin pixels (N s id c = 512). The power spec- 
trum was computed from a subset of the data on an 
clliptically shaped region of semi-axes a = 20 and 
b = 12 deg, which covers / s k y = 1.82% of the sky 
and is centered on the best observed region of the 
Boom-LDB flight. As shown in Fig. 1 the sky cover- 
age in this area, and therefore the noise per pixel, is 
nonuniform. The number of observations per pixel 
varies between 50 and 1510, with an average of 370, 
and a standard deviation of 120. The CMB dipolc 
was not included in the sky simulations as the high 
pass filtering used in map making reduces its rms 
residual variation in the observed region of the sky 
to less than 0.3 ^K, negligible compared to the rms 
amplitude of ~ 150 /iK of the intrisic small scale 
CMB fluctuations. 

The input CMB anisotropy spectrum used for the 
CMB sky model was chosen to fit the results of the 
joint Maxima-Boomerang data analysis (Jaffe et al. 
2001): fi = 1, A = 0.7, h = .82, n b h 2 = 0.03 and 



Apodisation W(r < p{4>)) w 2 /w4 
Top-hat 1 1 

Cosine cos ( 2 p(0) ) °- 514 
Gaussian exp^-jj^j^ 0.223 

Table 1: Apodisation applied to the clliptically 
shaped region of the sky used for extraction of 
the CMB anisotropy power spectrum. p{<fi) — 

a l \J 1 + (a 2 /b 2 — 1) sin 2 <j) is the radius of the ellipse 
in a direction <fi, and Wi — J dnW (n) 1 / An f s k y 

n s = 0.975. 

The number of MC simulations performed to es- 
timate the noise, the high-pass filter related trans- 
fer function, and the Ci statistics was, respectively, 
N$, = 450, N$ c = 250 and N^+ n) = 1350. Such 
high numbers of MC simulations are not necessary 
in practice, but were executed here to test accu- 
rately for possible biases, and to measure to a good 
precision the statistical distribution of the Ci esti- 
mates. Three different apodisations of the analysed 
sky region, described in Table 1, were used for power 
spectrum estimation. 

4.2. Statistics of Binned Spectrum Esti- 
mates 

Figure 5 shows the Monte Carlo averages and er- 
rors on the estimated C&-S, computed with the bin- 
width Al = 50 for the three used sky window apodi- 
sations. 

Clearly, the average recovered power spectrum is 
very close to the input one both at small £-s, where 
the signal is dominant, and at large is, where the 
signal to noise ratio is very low. This demonstrates 
that the estimator (26) is not biased. The error bars 
obtained from the Monte Carlo simulations (Eq. 29) 
can be compared to the "naive" ones (Eq. 16) cor- 
rected for the transfer function Fi. This effect de- 
creases the effective number of modes at each I to 

2 

^ = (24 + l)A^/sk y — F h , (35) 
which renders the new analytical error estimate 

V ' V V b 

Figure 5 shows that the MC error bars are almost 
identical to these analytical estimates. Figure 5 
also shows the spectrum estimate bin-bin correla- 
tion matrix, renormalised such that the diagonal is 
unity. This matrix is diagonally dominated and ex- 
cept for the first few bins, the off-diagonal elements 
are at most 10% as large as the diagonal elements. 
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Fig. 6. — Statistical distribution of the MASTER estimates of the binned full sky CMB anisotropy power 
spectrum are shown for the bins centered at I = 50, 200, and 750. The results of application of three different 
sky window apodisations, the Tophat, Cosine, and Gaussian, are shown in top to bottom panels, respectively. 
The abscissa shows the values of the binned spectrum estimates, C b , centered and normalised with the theoretical 
values used in the MASTER simulations: x = (C b -Cl h )/AC b (see Eq. 36). The histo gram shows the distribution 
derived from the MASTER simulations, the (blue) solid line is derived from the \t model where v is given by 
Eq. (35), and the (red) dashed line a Gaussian of same mean and variance. 



The distribution of the MASTER C b estima- 
tors measured in an £-b'm b defines the likelihood 
P(C b \C th ) of measuring C b given the input power 
spectrum C th (after likelihood marginalisation over 
all the other bins). This distribution was esti- 
mated from 1352 Monte Carlo simulations, and is 
illustrated in figure 6. The histograms show the 
simulated distribution together with the scaled xt 
model, where v (given by Eq. 35) is indicated on 
each panel, and a Gaussian distribution with the 
same mean and variance, which becomes undistin- 
guishablc from the \ 2 curve for large v. At all £-s, 
and specially for small values of i/, the \ 2 model, 
which has no free parameters, gives a much bet- 
ter description of the actual distribution than the 
Gaussian. 

To estimate the cosmological parameters one 
seeks a theoretical model Cf 1 defined by these pa- 
rameters, for which the likelihood P(C th \C b ) con- 



structed for a given data set is maximised. Using 
the Bayes theorem, this can be rewritten as 



P(C th \C b ) <x P(C b \C th )P(C 



■7th\ 



(37) 



where F(C th ) is the prior on the cosmological pa- 
rameters. The term P(C b \C th ) is often assumed 
to be a Gaussian function of the data C b , whereas 
its dependence on the theory C th can be approxi- 
mated by an offset-lognormal function (Bond, Jaffe 
& Knox 2000). We see however from Fig. 6 that 
for experiments with small sky coverage, the likeli- 
hood P(C b \C th ) has to be described at small is as 
\ 2 function of C b to avoid biasing the power spec- 
trum estimation and the cosmological parameters 
extracted from it. 



4.3. Monte Carlo Convergence of the Power 
Spectrum Estimator 

We checked, in the Monte Carlo simulations de- 
scribed above, that the estimators for the transfer 
function, the noise power spectrum (on the sky), 
and the accuracy of the error on Cb converge with 
the number of MC cycles, respectively, as follows 



sj\ 



a 10" 



100 



MC 



^ (38) 



with a between 0.6 and 0.9, depending on the sky 
window used; 



and 
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-44 ~ 0.6 10~ 2 
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U~\ T7> (39) 
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The contribution of the MC based estimation of F 
and M to the error on the recovered power spectrum 
is 

SCb _ 5Fi_ , 8M b ftb 
~~ ~ ~F\ 



C b 



+ 



Mb C b 



(41) 



and can be as low as a few percent for N$ c 



-^mc ~ ^00 in the signal dominated regime where 
Mb <C Cb- In the noise dominated regime, the ra- 
tio of the method induced uncertainty SCb to the 
statistical error ACb is 



SC b 
AC b 



0.L 



100 



(n) ' 
MC 



(42) 



can be brought below 10% in about 100 Monte Carlo 
simulations and analysis of pure noise TOD streams. 
Similarly, from Eq. (40), estimates of the indidual 
statistical errors of each binned ACb better than 
10% can be obtained with a few hundred MC cy- 
cles. We have seen however in section 4.2 that it 
is possible to predict analytically these errors with 
great accuracy. 



renders unbiased estimates of the CMB power spec- 
trum, with error bars very close to optimal. Monte 
Carlo calibration of the MASTER method requires 
generation and analysis of < 1000 independent sim- 
ulated realisations of a given CMB experiment. We 
demonstrated that CPU time requirements of the 
MASTER approach permit succesful analysis of the 
largest CMB data sets that exist at the present time 
using very modest computer facilities, for example 
an inexpensive PC farm. 

Because of the combination of the unsophisti- 
cated map making that we used and the aggressive 
high pass filtering applied to the Boomerang data 
stream in our numerical tests of MASTER, the es- 
timated power spectrum at the lowest multipolcs 
(£ < 100) has relatively large error bars. We are 
currently investigating the use of a more sophis- 
ticated map making algorithm to try to improve 
this situation. Ultimately, however, in the large sky 
coverage experiments, such as MAP or Planck, the 
power spectrum at low is can be analysed with fully 
hedged likelihood techniques if a coarsened pixelisa- 
tion is used. 

Possible improvements of the method include the 
modelisation of specific systematic effects, the use of 
more sophisticated map making techniques, and the 
extension to CMB polarisation measurements. 

EH would like to thank O. Dore for stimulat- 
ing discussions, J. Ruhl for useful comments on 
the manuscript, A. Lange for continuous encour- 
agement while this technique was developed, and 
all the Boomerang team for providing such a stim- 
ulating environnement. We thank A.J. Banday for 
help with creating the acronym for our method and 
for his careful reading of the manuscript. We ac- 
knowledge the use of HEALPix, cmbfast and fftw. 



5. Conclusions 

We have introduced a MASTER method for 
rapid estimation of the angular power spectrum of 
the CMB anisotropy from the modern CMB data 
sets. The method is based on direct spherical har- 
monic transform of the observed area of the sky and 
Monte Carlo simulations of the relevant details of 
observations and data processing. We demonstrated 
in an application of the MASTER method to simu- 
lated observations of Boom-LDB that this method 



A. Appendix: mode-mode coupling kernel 

In this appendix we compute the mode-mode coupling kernel resulting from the cut sky analysis, both in the 
planar geometry case, where the arithmetics involved may be more familiar and on the sphere. 

A.l. Analysis on the plane 

A scalar field AT(r) defined on the plane (or defined on the sphere and projected on a tangent plane) can 
be decomposed in Fourier coefficients as follows 

a(k) = J drAT(r)e- 2Mr , (Al) 

and 

AT(r) = / dka(k)e 2i7rkr . (A2) 



If AT is the homogeneous, isotropic, Gaussian distributed temperature fluctuation, each a(k) is an indepen- 
dent Gaussian random variable with 

(a(k)) = 0, (A3) 

and 

(o(k)o*(k / )) = 5(k - k')(C(k)) = <5(k - k')(C(k)), (A4) 

where 5 is the Dirac delta function. 

The Fourier coefficients derived on a weighted plane are then 

5(k) = J drAT(r)IT(r)e~ 247rkr , (A5) 

dk! a{k')K^[W]. (A6) 

If we write W(r) — J dkw(k)e 2z7rkr , the coupling kernel reads 

K klk2 = J dre 2l7rkir IT(r)e- 2l7rk2r (A7) 

= J k 3 w(k 3 ) J dre 2l7rr ( kl - k2+k3 ) 

= J k 3W (k 3 ) ( 5(k 1 - k 2 + k 3 ). (A8) 

If the polar coordinates of the vector k t are (ki,9i), the following useful property of the Dirac delta function 
can be demonstrated: 

J J d6 1 d6 2 5(k 1 + k 2 + k 3 ) = J d9 2 S(k 1 - |k 2 + kaD/fci 

= 27rJ(fc 1 ,fc 2 ,fc 3 ), (A9) 

where the function J is defined as follows 

J(k u fc 2 , h) = - {lk\k\ + 2k\k\ + 2k\kl -k\-k\- kj)~ 1/2 (AlO) 



for \k 2 — k 3 \ < k\ < k 2 + fc 3 , and J = otherwise. It follows that 

(■{k 2 + k 3 ) 
l(k 2 -k 3 ) 



/r(K 2 +K 3 j 
dk 1 J(k 1 ,k 2 ,k 3 ) = it I d(kf)J(ki,k 2 ,k 3 ) 
J(k 2 -k 3 ) 2 

= 2 y 1 du(i-u 2 y i/2 



= 2tt, (All) 

and J(k\, k 2 , 0) = S(ki — k 2 )/k\. 



(C kl ) = e»i<a(ki)o*(ki)), (A12) 

= ^tf de if rfk 2 / dk 3 (a(k 2 )a*(k 3 ))K klk2 [W}Kl ik3 [W] 
= ^f d0 if dk 2 {C k2 ) J d9 2 \K WlW2 [W}\ 2 



2tt J J k 2 dk 2 k 3 dk 3 (C k2 )W(k 3 )J(k u k 2 ,k 3 ), 



where W(k) = J d0w(k)w{k)* /2tt. 

This equation can be rewritten as follows 

(C kl ) = J k 2 dk 2 M klk2 (C k2 ), (A13) 

where the coupling kernel is given by 

M klk2 = 2tt J k 3 dk 3 W(k 3 ) J{k u k 2 , k 3 ). (A14) 

If C{k 2 ) = N (the white noise spectrum) then 

(C kl ) = 2irN J k 3 dk 3 W(k 3 ), (A15) 

where we used Eq. (All). 

A. 2. Analysis on the sphere 

A scalar field AT(n) defined on the sphere and weighted with an arbitrary window function W(n) can be 
expanded in spherical harmonics as follows 

a tm = J dnAT(n)W(n)Y e * m (n), (A16) 

= £>w [ dnY e , m ,(n)W(n)Y; m (n), (A17) 

t'm' 

= ^av m >Ki mVm ,[W], (A18) 

I'm' 

where the kernel K describes the mode-mode coupling resulting from the sky weighting. If W is z-axis az- 
imuthally symmetric, Ke m i> m > = Ke m i' m 5 mm ' ■ See Wandelt, Hivon & Gorski (1999) for an analytical calculation 
of Ki m i' m in the case when W is a tophat window. 

Note that a is a linear combination of Gaussian variables and is therefore Gaussian as well, but the de m s 
are not independent. 

If we use the series representaion of the window function, W(n) = w emYe m (n), the coupling kernel reads 
Ki imi e 2 m 2 = J dnY eimi (n)W(n)Y; 2m2 (n) (A19) 
= ^2we 3m3 J dnY timi {n)Y l:im:i {n)Yl 2 

7712 V 11 / 

"(24 + l)(2£ 2 + l)(24 + l) ll/2 



5>, 3 „ l3 (-ir 2 

^37713 



47T 



x , h e 2 e 3 \( ii ii h , (A20) 



where we introduced the Wigner 3-j symbol (or Clebsch-Gordan coefficient) y m m m y Several prop- 
erties of the 3-j symbol will prove useful. This scalar object describes the coupling of 3 angular momentum 
vectors (whose squared moduli are £i(£i + 1), and projections on the same axis are rm, for i = 1,2, 3) such that 

( h £2 h \ 

the total angular momentum vanishes. is non zero only if the triangle relation 

Y mi m 2 TO3 J 

\li-l2\<h<h+t2 (A21) 

is satisfied, and 

m 1 +m 2 + m 3 = 0. (A22) 
The orthogonality relations of the Wigner symbols read 



£3^.3 



(A23) 



T ( £l 12 h \( £l 12 ^=V* m! m'^i>^4)s7^-r, (A24) 

^ \ mi mo m,Q /I m,i m,o m/o / ^ m 3 m 3 L > - a ' ^ o/> i i ' v ; 



mi m 2 m 3 J \ m 1 m 2 mg / ^ "^"^ v ' ' '2£ 3 + \ 

where 6(£i,£ 2 ,£ 3 ) = 1 when the triangular relation (A21) is satisfied, and ^(^1,^2,^3) = otherwise. Finally, 
several recursive or closed form relations can be used to compute the 3 — j symbols. A useful example of the 
latter is 



h £2 t 3 \ = ( _ u L/2 

{ ' 



(L-2£ 1 y.(L-2£ 2 y.(L-2£ 3 )\ 



1/2 



(L/2)\ 



(L/2-h)\(L/2-e 2 )\(L/2-£ 3 )\ 



(A25) 



for even L = £\ + £ 2 + £ 3 (or equal to for odd L), with the asymptotic behaviour for L » 1 

(q 1 % ^ ^l(2£l£l + 2£l£i + 2£l£i-£\-£\-£f)- l/2 . (A26) 

See Edmonds (1957) for further details on Wigner symbols. 

The ensemble averaged power spectrum of the random scalar field AT(n) on the sphere computed with an 
arbitrary weighting function W(n) can be represented as follows 

1 ll 



1 11 

— E E E K™2 a lm 3 )^i m i^™2 [ H/ l if «im 1 l 3 m 3 ^ 
mi = —li £ 2 rri2 l^m^ 

— T E E^) E l^uWWll 2 - (A28) 

1 ' ™ n n n 
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mi =—ti t 2 

Upon substituting the kernel expansion in terms of Wigner symbols (A20), and reordering the sums, this 
expression expands to 

(C tl ) = E(^>^rE J2w e3m3 wl m4 ((2£ 3 + l)(2£ i + l)) 1/2 

£2 £3^3 £4^4 

£1 £2 £3 \ ( i\ £2 Mr/ *1 ^2 h\( £l £2 £4 



X 



/ V 1 \ mi —m 2 m 3 J \ mi —m 2 



y t! i 2 <* *i ^2 e 4 (A2Q) 
^ — ' \ m,i —mo m^ / \ m,i — m,o m,/i / 



which can be remarkably simplified with the aid of both the orthogonality relation of the Wigner symbols (A24) , 
and the definition (Eq. 10)of the power spectrum of the window function, We- The final expression reads 

(QJ = E%^), (A30) 



with 

"^^E^ + W? \ (A31) 

The Wigner symbols can be numerically computed from (A25) or from any equivalent recurrence relations. The 
equation (A30) expresses the ensemble averaged angular power spectrum measured with an arbitrary window on 
the sky for the statistically homogeneous and isotropic fluctuations described by an arbitrary ensemble averaged 
power spectrum over the full sky. 

If the input power spectrum is constant, (Ce) = N, corresponding to the white noise distributed over the 
sky, Eq. (A30) can be simplified using the orthogonality relation (A23), and the measured windowed power 
spectrum is also a constant 

(C e ) = Nf sky w 2 . (A32) 



B. Appendix: transfer function for parallel scans 

In the case of Boomerang, or of any scanning survey, a crude estimate of the transfer function corresponding 
to the high-pass filtering of the TOD, (see the Eq. 19), can be obtained by assuming that to first order 
the scans are parallel and performed at uniform angular speed. In such a case the filtering of the TOD will 
alter the sky signal structures parallel to scan direction, but leave unchanged the structures orthogonal to the 
scan direction. If the surveyed area is small enough (~ 20 deg in each direction) the tangent plane approach is 
sufficient to model the survey. Hence, the map can be decomposed in plane waves 

AT(x,y) = J a(k x ,k y )e^ x+k ^ (Bl) 

with k x — kcos9 and k y = ksuvO. The a(k) are zero-mean Gaussian variables with a variance 

(a(k)a(k , )*) = <P(fe))5(k-k') ) (B2) 
and the map power spectrum is given by 

p ( k ) = 77- I de a(k)a(k)\ (B3) 
2tt J 

If the scan is performed along the x-axis the map obtained from the filtered TOD is 

AT m (x,y) = J a(k x ,k y )f(k x )e^ x+k yy\ (B4) 

and its power spectrum is 

PfiitW = d6a(k)a(ky.f(k x ) 2 . (B5) 

The ensemble averaged filtered power spectrum is 

(Pmt(k)) = (P(k))^ J d0f(kcos9f. (B6) 

Hence, the effect of the TOD filtering on the power spectrum of the map amounts to a simple transfer function 
given by (£ <~ k) 

F £ (0) = ^- J d6f{kcos6) 2 . (B7) 

If the scan is performed at an azimutal speed v az at an elevation 9 e i , and the high pass filter applied to the 
data has a Gaussian form 

f(u) = 1 - e - {v/u " )2/2 (B8) 

then 

F^^2/tt r /2 ^(l-e-^ cos ^) 2 / 2 ) 2 , (B9) 



where i c — 2-KV c j(v az cos 9 e i). The asymptotic forms of arc 



*f - t«tc (BIO) 
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Ff - i-^J' <><c (BH) 
On the other hand, if the high pass filter is a chosen in the form of a sharp cut at the frequency v c 



f{u) = 0, v<v c , (B12) 

/M - 1, f>i>c, (B13) 

then 

Ff ] = 0, £<4, (B14) 

F £ (0) = 1-2/tt sin _1 (4/^), £>t c - (B15) 
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